COARSE-GRAINING MOLECULAR DYNAMICS MODELS USING 
AN EXTENDED GALERKIN PROJECTION 



XIANTAO LI * 

Abstract. We present a new framework for coarse-graining molecular dynamics models for 
crystalline solids. The reduction method is based on a Galerkin projection to a subspace, whose 
dimension is much smaller than that of the full atomistic model. The subspace is expanded by 
adding more coarse-grain variables near the interface between lattice defects and the surrounding 
regions. This effectively minimizes reflection of phonons at the interface. In this approach, there is 
no need to pre-compute the memory function in the generalized Langevin equations, a typical model 
of interface conditions. Moreover, the variational formulation preserves the stability of mechanical 
equilibria. 
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1. Introduction. Molecular dynamics (MD) models have been playing an in- 
creasingly important role in the study of micro- mechanical systems. The main dif- 
ficulty in the direct implementation of the MD model is the huge number of atomic 
degrees of freedom. Simulating a realistic system over a time scale of practical in- 
terest usually requires an extremely expensive computation. Therefore, it is of great 
interest to develop coarse-grained (CG) MD models in which only a small number of 
representative degrees of freedom are involved. For systems close to equilibrium, a 
coarse-graining energy can be obtained by integrating out the remaining degrees of 
freedom in the free energy. This is the idea behind most direct CG models, such as the 
one by Rudd and Broughton [64, 63] for crystalline solids. The equations of motion 
is then expressed using the Hamilton's principle. Notice, however, these models are 
not derived directly from the MD model. In particular, for systems far from equi- 
librium, e.g.^ propagation and interaction of shock waves, this procedure can not be 
applied. There have also been attempts to derive CG models using homogenization 
methods [17, 24, 34]. Such asymptotic approach offers a means to derive a continuum 
mechanics model, although such procedure can only be applied to perfect crystals. 
In contrast, Seleson et al [66] considered the upscaling to peridynamics (PD) models 
[67]. 

Another approach that has gained great popularity is to directly combine MD with 
continuum mechanics models. The idea is to retain the atomic-level description near 
critical areas where the molecular trajectories are of interest, while in the surrounding 
region, the MD model is replaced by a continuum mechanics model, often in the form 
of partial differential equations. At the interface, a coupling condition is imposed 
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to combine the two models, see [1, 27, 28, 62, 18, 74, 77, 61, 72, 29, 50, 69, 70, 59, 
76, 60, 8, 52, 43, 33]. Despite the many proposed methods, quite a few important 
issues have not been fuhy addressed. Examples include systems at finite temperature, 
dislocation dynamics problems where dislocations pass through the interface, and 
adaptive selection of the MD region. In addition, stability and error analysis for such 
coupling methods is open. 

In this paper, we will consider the Galerkin projection procedure to derive CG 
models for crystalline systems. The method is motivated by the Galerkin approx- 
imations of some Hamiltonian PDEs, especially the second-order wave equations. 
Examples include the finite element approximations of the linear second-order hyper- 
bolic equation [10, 11, 26, 35] and nonlinear elasto-dynamics equation [54]. In fact, 
as demonstrated in [5, 30, 66], for perfect crystals the MD model at zero temperature 
can be approximated by the elasto-dynamics models with the constitutive relation 
given by the Cauchy-Born rule. Within the conventional Galerkin framework, the 
coarse-grained MD model can be derived by projecting the equations of motion to the 
subspace spanned by the basis functions. Near defects, this projection leaves the MD 
equation unchanged, while in the surrounding region, the CG equations coincide with 
those in [75], obtained by a least square procedure. The CG model obtained from this 
conventional Galerkin projection is also similar to the model developed by Chen and 
her coworkers [18, 19, 78], which were derived from the Irving-Kirkwood formalism 
[39]. But in their model, the issue of phonon refiection has not been addressed. 

As we will show, the conventional Galerkin method often leads to significant 
phonon refiections at the interface. The refiection can be attributed to the discrepancy 
between the dispersion relations in the atomistic and CG regions, as well as the varying 
mesh size. In fact, when solving wave equations, non- uniform mesh often causes inter- 
grid refiections [14, 15]. 

To minimize the phonon refiection, Wagner, Liu and their coworkers [75, 69, 70] 
developed a bridging scale method in which the system is divided into a coarse scale 
region and a MD region. In the coarse scale region, the displacement of the atoms 
consists of a coarse scale part, the evolution of which is governed by equations de- 
rived from the Galerkin projection, and a fine scale part, which can be eliminated by 
linearizing the MD model and using Laplace transform. This leads to a coupling condi- 
tion at the interface, expressed in the form of a generalized Langevin equation (CLE). 
The reduction procedure is closely related to the derivation of absorbing boundary 
conditions for MD models [3, 2, 73, 74, 41, 42, 51, 47]. In order to implement this 
coupling condition, the memory kernel function in the generalized Langevin equation 
has to be precomputed. This calculation is a highly nontrivial task, especially because 
the functions tend to have very slow decay and they depend on the geometry of the 
interface. This numerical issue has been addressed in [27, 28, 52], where the memory 
function is approximated by functions with compact support. A similar approach 
to the BS method is the perfectly matched (PMM) multiscale method [72, 46], in 
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which the MD equations are modified by adding damping terms and modifying the 
force constant matrices at the interface to prevent artificial refiection. This method 
is easier to implement. However, it is not clear how this can be correctly formulated 
in the case of coupling methods. Ideally, one should only remove the waves modes 
that are not captured by the coarse scale variables. But no such selective damping 
mechanism is present in the design of the PMM method. 

A common open problem in the methods mentioned above is the well-posedness of 
the coupled models. By imposing additional conditions at the interface, the variational 
structure is lost. In fact, for an isolated atomistic system, the stability of the boundary 
conditions is already nontrivial [48]. We will show some examples in this paper that 
some truncation schemes may lead to unstable models. In this paper, we propose a 
new approach to extend the conventional Galerkin method. The idea is to expand the 
approximation space while maintaining the variational formulation. We start with a 
conventional Galerkin method. This is followed by a few refinement steps. At each 
step, we keep the mesh structure intact. Rather, we define additional variables at each 
nodal point. Therefore there is no need to re- mesh frequently to enhance the overall 
accuracy. The coarse-grained MD model is derived by projecting the equations of 
motion to the extended subspace. Therefore the variational structure is kept. As the 
expansion continues, we are building a hierarchy of ODE systems, which in the limit 
is equivalent to the full atomistic model provided that the initial set of basis functions 
are chosen properly. In particular, the conventional Galerkin approximation is at the 
top of the hierarchy. It can be obtained by dropping all the additional variables. From 
this perspective, the large error that usually arises at the interface can be attributed to 
the premature truncation of the hierarchal structure. This structure is reminiscent of 
the BBGKY approximations of the many-body Liouville equation in non-equilibrium 
statistical physics [12, 16]. Our expansion procedure is in the spirit of higher moment 
models, e.g.^ the Bennett model [4]. By introducing such higher order models at the 
interface, we are able to minimize phonon refiections and improve the accuracy of the 
CG model. 

This Galerkin projection approach bears similarity to the Mori-Zwanzig (MZ) 
projection formalism, a well known reduction methodology in nonlinear statistical 
physics [57, 56, 82]. For systems out-of-equilibrium. The MZ formalism has been 
used to derive coarse-grained molecular dynamics models [58, 25, 40, 49], and it has 
proven to be a very useful way of thinking about this type of problems. In the MZ 
formalism, one uses a projection operator to separate out coarse-grain variables, and 
then derive the effective dynamics. In the works on Zwanzig [82] and the recent work 
of Chorin and his coworkers [21, 20, 23, 22], conditional expectation is used as the 
projection. This is convenient when the initial data are sampled from an equilibrium 
distribution, e.^., the canonical ensemble. In principle, one can use this approach 
to derive a GLE as a CG model [49]. But it suffers the same problem mentioned 
above - the computation of the memory kernel function. In fact, the expression of 
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the memory function involves the inversion of a large matrix, with dimension equal 
to the total number of degrees of freedom in the system. On the other hand, crude 
approximations, e.g.^ by Dirac delta functions, are too ad hoc, and often introduce 
large modeling error. 

The rest of the paper is organized as follows. We first describe the conventional 
Galerkin method, applied to the MD model. We then show a simple numerical test 
to demonstrate the accuracy of this method. Section 3 presents the main idea of 
extended Galerkin method, followed by some tips for practical implementation. In 
section 5 we give some results of the numerical experiments. 

2. The Conventional Galerkin Projection Method. 

2.1. The mathematical formulation. The formulation has strong analogy to 
the Galerkin projection of second order wave equations. But there are also crucial 
differences that will be emphasized. We consider a system of atoms with reference 
position in domain Q with boundary dQ. A schematic illustration is shown in Fig. 
2.1. 







• 


• 








o 


o 






• 


• 


• 


o 


o 


o 


o o 










o 


o 


o o 


o 


o 


o 




o 


o 




o 


o 


o 




o 


o 


o o 


o 


o 


o 




. o 


o 


o o 


o 


o 


o 


• 




o 


o o 


o 


o 


o 



Fig. 2.1. A schematic illustration of the domain. Open circles: the atoms in the interiors- 
Closed circles: the atoms at the boundary. 

We let Yi and be respectively the reference and current position of the ith atom, 
with its displacement given by = — y^. For the molecular dynamics model, the 
atoms obey the Newton's second law, 

m^Xi = -Vx.K (2.1) 

Here rrii is the mass of the ith atom, and x^ denotes the second derivative in time. 
Further V = y(xi,X2, • • • ,X7v) is the interatomic potential, which will be assumed 
to be an empirical model. The boundary condition is assumed to be of Dirichlet type, 
i.e., the position of the atoms in dft is prescribed. 
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We first write the equations of motion in the following compact form, 



u = f(u), f = -VV. (2.2) 

Here we have chosen to work with the displacement and grouped the displacement of 
all the atoms into one vector. We have set the mass to unity and neglected possible 
external forces for simplicity as they will not significantly change the formulation. 

To explain the numerical procedure, it is often helpful to separate out the har- 
monic and anharmonic parts of the interaction as follows, 

f(u) = -Au + fAr(u). (2.3) 

Here ^N indicates the anharmonic part of the force. The harmonic interaction can 
often be defined by a linearization around the reference configuration. See [79] for 
such examples. We will refer to the model 

ii = -Aw, (2.4) 

as the linearized MD model. 

Next, let n be the dimension of u, and let X = be the entire space. To con- 
struct a space for the approximate solutions, we choose a set of m independent nodal 
basis vectors, < i < m}. Let ^ = (fi, (p2^ ' - - ^ be a matrix, each column 

of which is a basis vector. We denote Y = Range($) the m-dimensional subspace. 
The matrix ^ may be constructed from an interpolation procedure for continuous 
problems, and each function : ^7 ^ R is a shape function. One can also think of 
the matrix $ as an interpolation operator from R^ to X, and similarly, can be 
considered as a restriction operator, similar to the terminology used in algebraic multi 
grid methods [55]. To better illustrate this idea, we give two examples. The first ex- 
ample is similar to a finite element representation with piecewise linear functions. The 
top figure in Fig 2.1 shows one such function in the subspace. In this case, the nodal 
points coincide with the reference position of some atoms. The middle figure shows a 
similar function, but the nodal points are away from atom positions. For the second 
example, we consider a piecewise constant function, which is similar to first order 
finite volume representation. We remark that the basis functions are discrete func- 
tions, although they can often be obtained by confining a continuous basis function 
to atomic positions. Near lattice defects, we often retain all the atoms in a neighbor- 
hood. For these atoms, we choose the trivial basis function which equals to one at an 
atomic position and zero at all other atoms. This part of the construction is due to 
the discreteness of the MD model and has no analogy in Galerkin approximations of 
PDEs. 

In X, we define the usual inner product, 

(u,v) =^u,.v,. (2.5) 
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Fig. 2.2. Interpolation of a discrete function. Top: A piecewise linear interpolation in which 
the nodal points coincide with the reference position of some of the atoms; Middle: Another piecewise 
linear interpolation where some of the nodal points are away from the atomic position; Bottom: A 
piecewise constant interpolation. 

In the Hilbert space, we define the orthogonal projection to the subspace which in 
matrix notations, can be expressed as 

P = ^{^^^)-^^^ . (2.6) 

In addition, we let Q = / — P be the complementary projection. 
In addition to the inner product, we define, 

a(u,v) = -(f(u),v). (2.7) 

For linear models, this is reduced to, 

a(u,v) = n^Av. (2.8) 

Using some symmetry properties of A, this quantity can often be written as [7], 

a(u, v) = ^(u^ - WjYAij{Yi - Vj). 

ij 

Where Aij is a force constant matrix. This quadratic form is similar to the bilinear 
form of elliptic PDEs. 

We first discuss the conventional Galerkin approximation. Such an approach seeks 
an approximation u : [0, T] F, such that, u G i^^(0, T; X), and for any test function 
V G F that vanishes at the boundary, there holds, 

(G,v) +a(G,v) = 0. (2.9) 
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If we let M — q = M ^$^u, then u = ^q, and we have, 

Mq = <l>^f(<l>q). (2.10) 
For the Unear problem, these equations become 

Mq=-Kq, (2.11) 

where 

K = ^^A^. (2.12) 

The matrices M and K can be regarded as the mass and stiffness matrices in finite 
element approximations of elasto-dynamics problems [37]. 

2.2. An error analysis. The induced norm will be denoted by || • ||o. For a 
given time interval [0,T], we also define. 



|u||l-(0,T;X) = sup ||u(t)||o, 

te[o,T] 



and. 



|u||l2(0,T;X) 



Jo 



In addition, we define the Ritz projection operator to the subspace, 

P = ^{^^A^y^^^A, (2.13) 

with Q = I — P being the complementary projection. 

To compute the error, we let e = u — u = 6> + ?7, in which. 



= u - Pu, 
77 = Pu — u = 



(2.14) 



For the initial condition, we will assume that 

u(0) = Pu(0), ii(0) = Pu(0). (2.15) 

The analysis is similar to the analysis of finite element methods for second-order 
wave type of equations [10, 11, 26, 35, 45]. The analysis presented here follows the 
approach in [10], and provides the error estimate in the || • ||o norm. 

We consider the Galerkin method for the linear problem. In this case, the error 
equation for in the variational form is given by, 

{0,y,)+a{0,y,) = -{i^,^). (2.16) 



for any cp gY. 
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We now re- write the equation as, 

^ ((9, cp) - ((9, cp) = -a{0, ^ ^) + ^) • 
Next fix a r > 0, we choose a particular test function as, 

ip{t) = f 0{s)ds, (2.17) 



Jt 

This special choice implies that 

(fir) = 0, ip = -e. (2.18) 
With this choice the equation above becomes, 

|i(M) = -|(e.^) + |l«(^,„) -(«,(.). ,2,<,) 

Integrating the equation from to r, one gets, 

= \\emi-a{m,m) 

-2(e(0),^(0)) -2 r {fj{s),e{s))ds. ^^'^'^^ 

^0 

We will choose the initial condition so that (e(0), (/:?(0)) = 0. In addition, we notice 
that, a((/:?(0), ^(fS)) > 0. Therefore, using Cauchy- Schwartz inequality, we get for each 
r G [0,T], 

mr)\\i<\\oml^\\mL-io,T^^^^ (2.21) 

This leads to, 

m\L^iO,T;X) < \^ll^(0)||o + 2Vf\\ri\\L2^o,T;X)- (2.22) 

Combining with (2.14), we get. 

Theorem 1. The error of the Galerkin method for the linear problem satisfies, 

||U - u||^oo(o,T;X) < C'(||gu||^oo(o,T;X) + VT||Qu||l2(0,T;X)) • (2.23) 

Remark 2. For perfect crystals, it is often possible to introduce norms that are 
similar to H^{Q) and H'^{Q) norms. In this case, one might be able to derive a 
more explicit bound for the interpolation [53] and obtains an error estimate in H^{Q) 
norm. In general, however, some symmetry may be lost near a lattice defect, it would 
be difficulty to define such norms. 

Remark 3. There are two cases when such an estimate loses its values. The 
first case is when the time scale is beyond C^(l) scale, i.e., the molecular times scale, 
which is about the scale of pico- seconds (l{)~^'^s). In fact, long time simulations of 
molecular systems is needed for most problems. Another problem is when phonons are 
generated by local defects. Due to high frequency phonons, the term \\Q\^\\ might be 



2.3. A one-dimensional example. To illustrate the idea of using the Galerkin 
projection as a reduction method, let's consider a one-dimensional model. 

mxj = V\xj^i - Xj) - V'{xj - Xj-i). (2.24) 

We use the Lennard- Jones model for the potential energy V, given by, 

V{r) =r-^^ -r-^ (2.25) 

The lattice spacing ao and the mass are both normalized to unity. We consider a chain 
with 1024 atoms. We will coarse-grain the first 512 atoms and keep the remaining 512 
atoms. This divides the system into two sub-domains. In the coarse-grained region, 
we choose one node out of every eight atoms, and we use usual piecewise linear basis 
function. On the other hand, in the MD region, every atom is chosen as a node. We 
impose homogeneous boundary conditions at the ends. For the initial condition, we 
start with a left-moving wave packet centered in the domain on the right. To create 
such initial condition, we first define, 

20 20 

u\x) = ^cos(a^), v\x) = ^-c^(a)sin(a:r), (2.26) 

where uj{C) — 2 sin | is the dispersion relation, and = 0.5 + 0.02/c selects a few wave 
numbers. Then we choose the initial condition as follows, 

u,(0) = F(j>0(j), t;,(0) = FO>0(j). (2.27) 

Here F{x) = 0.00025 exp{ — (x — 640)^/800} is a Gaussian profile to confine the wave 
packet inside the atomistic region. 

The top figure in Figure 2.3 shows the result from a full atomistic simulation. 
We observe that the wave packet moves across the interface and enters the region on 
the left. For the conventional Galerkin method, however, the wave packet arrives at 
the interface, and then most of it is reflected back into the region on the right. We 
may attribute the reflection to the large ratio between the mesh size and the atomic 
spacing. Therefore, we conducted another test using the Galerkin method in which 
the mesh size changes linearly from 8ao to ao as we move from the interface to the 
left boundary. The figures on the third row in Figure 2.3 show the corresponding 
results. We see that although the mesh size varies very slowly, we still get significant 
reflections. Compared the case with uniform mesh size, the reflection occurred a bit 
later. 

3. The Extended Galerkin Method. We have observed from the one-dimensional 
example that the conventional Galerkin method may give poor accuracy, especially 
when phonons are propagating out the atomistic region. In this section, we aim to 
improve the CG model (2.10) derived from the conventional Galerkin method. 
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Fig. 2.3. Numerical results based on the one- dimensional model. 



Fig. 2.4. A numerical test with the one- dimensional model. From top to bottom: The results 
from a full atomistic simulation; The results from the conventional Galerkin method; The results 
from the conventional Galerkin method with slowly varying mesh size. From left to right: Snapshots 
at T=20, 40, 60 and 80. 

3.1. Galerkin approximation with memory. We first derive an exac^ model 
for the coarse-grained variable q. Let y{t) = ^^\i(t). We set our goal to finding 
an effective equation for y. To better illustrate the ideas, we first consider the lin- 
ear model. Taking time derivatives and split the right hand side according to the 
projection operators, we get. 



y = - 



(3.1) 



The first term on the right hand side can be written as. 



Here M and K are the mass and stiffness matrices defined in (2.11) and (2.12). 
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Notice that by dropping the second term on the right hand side, and letting 
q = M~^y, we obtain the conventional Galerkin approximation (2.9) or (2.10). The 
example from the previous section suggests that a truncation at this step would be 
premature. To obtain a more accurate solution in time, we do not completely neglect 
the second term. Instead, we define, 

z = AQu. (3.2) 

Again taking derivatives and split the right hand side in a similar way, we find, 

z = -AQz - AQAPu = AQz - AQA^M'^y. (3.3) 

To explicitly express z, we define fundamental solutions, C{t) and S{t)^ which satisfy 
the equations respectively, 

C = -AQC, C{0) = /, C(0) = 0, (3.4) 

and 

S = -AQS, S{0) = 0, S{0) = I. (3.5) 

Clearly, these two functions both commute with AQ. In addition, they are related as 
follows, 

S{t) = C{t), C{t) = -S{t)AQ. (3.6) 

Formally, the two functions can be written as, 

C{t) = cos(AQ) + t, S{t) = (AQ)-+ sin(AQ)it. 

Although the matrix AQ may not be invert ible, the second expression can be defined 
by the corresponding power series of the sine function. 
With these two functions, we get, 

z(t) = C(t)z(0) + 5(t)z(0)- / S{t-T)AQA^M-^y{r)dr. (3.7) 

Jo 

By collecting terms, we find that, 

y = -$^APu+/ ^^S{t-r)AQA^M-^y{r)dr^^^g{t), (3.8) 
Jo 

where, 

g{t) = - [c(t)AQu(0) + 5(t) AQu(O)] . (3.9) 
Finally, we let q = M~^y and we obtain, 

Mq{t) = -Kq{t) ^ [ 6(t - r)q(r)dr + ^^g(t). 
11 



(3.10) 



We let w = Pu = $q. It satisfies an equation in a weak form, 

(w,v)+a(w,v)= / (P(t-r)w(r),v(t))dr+(g,v). (3.11) 
Jo 

Here the memory kernel functions are given by, 

e{t) = ^^B{t)^, B(t) = S{t)AQA. (3.12) 

O is a symmetric, matrix- valued function. 

The equation (3.10) is an exact equation for q. We now make an approximation 
by neglecting g(t). To justify this step, and to find more explicit expressions of the 
kernel function O, we state the following lemma. 

Lemma 4. Let A G cr{AQ) be an eigenvalue of AQ, and ^ he the corresponding 
eigenvector. Then 

1. A = O^^Gr. 

2. If X ^ 0; then X is an eigenvalue of QAQ, and the corresponding eigenvector 
is T] = Q^. 

3. If X is an eigenvalue of QAQ, and the corresponding eigenvector is 
T] G then X is an eigenvalue of AQ, and (, = r] -\- X~^PAr] = X~^Ar] is 
the corresponding eigenvector. 

4. Let Xi^i = 1, 2, • • • , be the eigenvalues of QAQ, and rji he the corresponding 
eigenvectors with unit length, then 

QAQ = ^Xir],r]f, Q = Yl = Yl ^^^^^^ 



5. The force g{t) is expressed as, 

g{t) = ^^A Y cos ^/Xit TjiTjf u{0)^^^ A Y 



sin vA7t 



A,Ga(QAQ)\{0} A, Ga(QAQ)\{0} 



T 



(3.13) 

6. The history function Q{t) can be expressed as: 

G{t) = ^^S{t)AQA^ = Y '^^^^{^^Arj,){^^Arj,f. 



XieaiQAQ)\{0} ^ ' 

Notice that the zero eigen-modes of AQ do not appear in the above expressions 
because of the orthogonality conditions. 

As a result, when the initial condition is close to the subspace F, g(t) would be 
small. Therefore, we may neglect g(t) to obtain a closed system. In particular, when 
u(0) ± Y and u(0) ± F, g(t) = 0. 

This yields the following equation, 

Mq{t) = -Kq{t)^ [ e{t - r)q{r)dr. (3.14) 
Jo 
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The solution u is approximated by u = $q. The corresponding weak form is given 

by, 



(ii,v)+a(u,v)= / (e(i-T)u(T),v(t))dT, 
Jo 



(3.15) 



for each v(t) G Y with homogeneous boundary condition. For nonhnear problems, 
one may replace the bilinear form a(u, v) by the nonlinear form (2.7) and express the 
effective equation in a similar form as (3.14). 

We have computed the memory function O for the one-dimensional example de- 
scribed in the previous section. Figure 3.1 shows the diagonal entries of the function 
Q{t) for points near the interface. We can make the following observations: The 
kernel function exhibits a peak near the interface, and away from the interface, the 
kernel function becomes significantly smaller. In particular, inside the MD region, it 
becomes zero. 
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n=62 



n=63 



n=64 



n=65 



Fig. 3.1. The diagonal entries of the kernel function, On,n{t), for the one- dimensional example, 
with n = 62, 63, • • • , 67. The bottom figure shows the corresponding nodes, n = 65 is the location of 
the interface. 



3.2. Approximations of the GLE with direct expansion. The GLE model 
provides a more accurate approximation by including a history-dependent term. In 
practice, however, the memory function has to be obtained prior to the computation. 
In general, direct calculation is quite difficult since the explicit expression involves 
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solving (3.6) which is almost as complicated as solving the original problem. Therefore 
the GLE often has to be approximated. One approach is to extend the system by 
introducing extended variables and embed the GLE into a system with no memory. 
This has been done for colored noise systems where the power spectrum is a rational 
function [13]. For our problem, we may rewrite the equation (3.1) into 

y = -KM-^y - ^^AQu. (3.16) 

The second term on the right hand side is not in the subspace Y. To proceed, we 
define, 

Co = ^^AQu. (3.17) 

Taking the second time-derivatives, and rearrange the terms on the right hand side, 
we can obtain, 

Co = - [K2 - K^K^^K,]q + i^ii^o"'Co - ^^A^Qu. (3.18) 

Here, we have defined generalized stiffness matrices, 

= ^^A^^, (3.19) 

for any ^ > 0. Notice that Kq = M, and its invertibility follows from the independence 
of the basis functions. 

This equation contains a term that can be represented by neither q nor Co- So 
we define another additional variable, 

Ci = ^^A^Qu, (3.20) 

This new variable satisfy the equation, 

= - [Ks - i^2i^o"'^i]q + K2M^Ho - KiM^Hi - ^"^AQA^Qn. (3.21) 

The last term in this equation is a new quantity. Conceivably, this process can be 
continued indefinitely, and a hierarchy can be built. This is similar to the BBGKY 
moment system for the Liouville equation. To obtain a closed finite system, truncation 
is often needed. For instance, we may drop the last term in the equation above, and 
we find an extended system, 

i^oq = -^iq-Co. 

< Co = -[i^2-i^ii^o"'^i]q + ^i^o"'Co-Ci, (3.22) 

We write this system in a more compact form, 

Mq = -^q, (3.23) 
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where, 





q 




Ko 








q = 


Co 


, M = 





/ 







. . 










/ 



K = 



Ki 10 



We further remark that, 

Remark 5. Such extension is not unique. For instance, one can instead define 
the additional variables as follows, 

Co = AQu, = {AQfu, ^2 = {AQfu, 

and the extended system reads, 

io = -(^2-i^ii^o"'^i)y-Ci, 



Remark 6. The additional variables may not have direct physical meanings as 
the system is further extended; 

Remark 7. The variational structure is lost, i.e., there is no weak form associated 
with the truncated system, and the extended system may not be well-posed. 

To demonstrate the potential ih-posedness of the extended system, let's consider 
the one-dimensional chain model with nearest neighbors. We choose a chain of 5 atoms 
with indices from to 4, and we fix the first and the last atom to their equilibrium 
position, i.e., i^o = 0,tt4 = 0. With some normalization, we choose the matrix A as 
follows. 



2 


-1 








-1 


2 


-1 








-1 


2 


-1 








-1 


2 



This matrix is positive definite. Therefore the trivial solution of the linearized MD 
model is stable. We now choose one CG variable with nodal value defined at the third 
atom, and we set, 

which corresponds to a linear interpolation. The stiffness matrices can be easily 
computed. They are given by, 

= Ki = K2 = 1. 

We consider the extension to by dropping from the equation (3.18). In this 
case, we found that the generalized eigenvalues for {K^ M) to be ±a/6/3. The negative 
eigenvalue indicates instability. 
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3.3. An extended Galerkin's approximation. From the example at the end 
of last section, one finds that such direct expansion is not reliable since it may in- 
troduce unstable models. We propose a new extension method which maintain the 
Galerkin formulation. This extension is done by expanding the approximation space. 
This would eliminate the problem of instability. One idea would be to extending 
the subspace by introducing bubble elements [6, 9]. In particular, this idea has been 
applied to multi scale problems [38]. But here we will follow a different path. 

For the purpose of the extension, we let Yq = Range(^). We let Yi = Range(QA$), 
and Y = Yq ^Yi. Then the Galerkin projection of the linear problem to Y yields 




-^2q-^3Co. 



(3.24) 



Here, the approximate solution is expressed as, u = + QA^^q. The matrices K2 
and Ks are given by, 

K2 = ^^AQA^ = K2- KiK^^Ki, (3.25) 
Ks = ^^AQAQA^ = Ks- K2K^^Ki - KiK^^K2 + KiK^^KiK^^Ki. (3.26) 

We now go further to extend the approximation space. For example, we can 
introduce another projection. 

Pi = QA^{^^AQA^)-^^^AQ, Qi = I - Pi. (3.27) 

We let Y2 = Range(Q(5iA^<l>). One can easily verify that Y2 ± Yq and Y2 ±Yi. We 
set F = Yq ® ^1 ® ^2 , and seek the solution in the form of, 

u = ^q{t) + QA^^oit) + QQiAH^^it). (3.28) 

Using the variational approach, we find. 



Koq = 


-K,q- 






mo = 


-K2(l- 




(3.29) 


. ^4^1 = 









Here the matrices are given by, 

K4 = ^^AQAQQiAQA^, 

^5 = ^^AQAQiQAQQiAQA^. 

For nonlinear problems, the Galerkin projection will be applied to (2.2). The varia- 
tional procedure yields, 
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For the linearized model, the coarse-grained model derived from this extended 
Galerkin method can be viewed as an approximation of the generalized Langevin 
model (3.14). To see this, we assume that u(0) G Y and u(0) G F, and in light of 
this, we set ^q{0) = 0. Applying Laplace transform to the second equation in (3.24), 
we found that, 

Jo 

which can be substitute into the first equation, yielding, 

Mci=-Kq^[ Slit- s)q{s)ds. (3.30) 
Here the function 0i has Laplace transform, 

ei(A) =K2(A'^2+^3)"'^2. (3.31) 

Similarly, (3.29) can be turned into the same form with the kernel function's 
Laplace transform given by, 

62(5) = K2 \{s^K2 + Ks) - {s^K^ + ^5) "'^4] K2. (3.32) 

This provides a rational approximation of the Laplace transform of the kernel function. 

We briefly explain the motivation for using such subspaces for the extension pro- 
cedure. Notice that when the initial data are in the subspace F, the solution of 
the harmonic model (2.4), is in the space, K = span{A^<l>, £ > 0}. This can be 
seen by writing down the solutions in the form of sine and cosine functions, and 
then expanding the solution into power series. Therefore, we may consider the sub- 
space Y = Range($) + Range(A$) = Fq © ^i- Similarly, one can directly verify that 
Y = Range(<l>)+Range(74^)+Range(74^<l>) = Yo^Yi^Y2. To continue with the expan- 
sion, one may consider the sequence of subspaces Kl{A^ ^) — span{<l>, A^^ • • • , A^<l>}, 
which is a Krylov type of subspace. The Cay ley-Hamilton theory implies that C 
Kn for ^ > n; n is the dimension of the problem. If Yq is not contained in the sum 
of the eigenspaces associated with a subset of the eigenvalues, then = X, and the 
exact solution should be recovered if the subspace Y is expanded to Kn. 

3.4. The one-dimensional chain problem. We now turn again to the one- 
dimensional example. As in the conventional Galerkin method, we partition the 
system into a CG region where the mesh size H = 8ao. For the extended Galerkin 
method, we need to select additional nodes near the interface. We choose 20 nodal 
points in the CG region, and 5 atoms from the atomistic region. The results are 
included in Fig. 3.4. For the choice L = 5, some reflections are observed. When L is 
increase to 10, almost no reflection is observed. 
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Fig. 3.2. The solution of the Id chain model using the extended Galerkin method. Top: L = 5; 
Bottom: L = 10. 

4. Implementation Issues. In the formulation of the extended Galerkin method, 
we have assumed that all the quantities, e.g.^ the projection matrices and the stiffness 
matrices, can be computed exactly. Clearly, computing these quantities exactly would 
require an effort that is as much as solving the original problem. 

4.1. Construction of the approximation spaces. We first discuss the con- 
struction of the subspace to which the Galerkin method is applied. Several steps are 
taken to simplify the procedure. 

First we discuss the subspace Yq. In light of the lattice defects in the system, 
where the molecular dynamics model should be kept, it is convenient to partition the 
system into a coarse-grained region, and an atomistic region. In the CG region, we use 
a conventional finite element representation of the displacement field, e.g.^ piecewise 
linear functions. Inside the atomistic region, we pick the trivial basis function, (p{'x.) = 
6nk at the point y^. This function is one at this point, and zero at all other points. 
The union of such point will be donated by l^at- From the previous section, we see 
that the memory function Q{t) is small away from the interface. Therefore, we will 
only select the basis functions with nodal points near the interface to expand the 
subspace. For each selected basis function, we let 

Kl{A, cp) = span {cp, Acp,-- - , A^cp} . (4.1) 

The whole extended subspace is constructed by combing K l{A] ip j) j G J, where J 
is the set of selected basis functions. For this purpose. We need to include some more 
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atoms around l^at- The new set is denoted by l^at- This wih be explained next. 

We now discuss how to generate Kl{A^ cp). We need to compute a vector in the 
form of AiIj where ip is the discrete function whose support hes in l^at- The vector 
win be approximated by a finite difference formula, 

f(y + e^)-f(y) ^ 

e 

Here y is the reference position. For most empirical potentials, the interaction has 
finite range, and the shape function has support in l^at, also has support close 
to Atp. In our implementation, we will choose l^at to contain the support of all 
A^(Pj,£ < LJ e J. 

Finally, the new set of functions generated from the previous step may not be 
independent. One implication is that the mass matrix may be singular. In this 
case, we will use an orthogonalization procedure similar to the Anoldi's method [65]. 
However, we only make the additional basis function orthogonal to Yq. This way, 
only the position of the atoms in (^at are needed in the process. We will denote the 
additional set of basis functions as, 

^ = [V^i,V^2,--- (4.2) 

This is a set of orthonormal basis vectors. In addition, "^^^ = 0. With this notation, 
we have F = Yq © Range(^), and the Galerkin projection yields, 

Mq = $^f(u), 

where u = + 

4.2. Summation and quadrature rules. The coarse-grained MD model is an 
autonomous, second order ODEs. We use the Verlet's method to integrate the system 
in time. One advantage of this is that at each step of the time integration, the forces 
only have to be sampled once. 

We now discuss the calculation of the mass matrix and the averaged forces. As 
an example, we consider the case where the CG region has been divided in triangles 
(2D) or tetrahedrons (3D), with basis functions cpn given by standard piecewise linear 
functions. For the mass matrix, we need to compute 

Mm,n = ^rn{yi)^n{yi)^ (4.4) 
i 

and for the forces, we need to compute, 

F„ = ((p„(y),f(y)), (4.5) 

where cpn is one of the basis functions. 

For the nodal functions in <l>, they are divided into two groups, the ones in ^lat, 
which includes the basis functions in ^, and the ones in the coarse-grain region. For 

19 



the formulas (4.4) and (4.5), if the basis functions have support in l^at, we wih evaluate 
the expression exactly. Otherwise, we will approximate the summation by integrals. 
This will be discussed next. 

We first consider the computation of the mass matrix. Let m and n be two vertices 
of a triangle T (or a tetrahedron in 3D). The sum in (4.4) is over the intersection of 
the support of cpm and cpn. It suffices to consider the sum over a triangle, which can 
be approximated by an integral, as was done in [75], 



Here Vo is a volume of the unit cell. This integral then can be further approximated 
by standard quadrature formulas. 

Next we consider the sampling of forces, given by the inner product. 



Quadrature approximation of such average has been discussed in [75, 36]. Here we 
will derive an expression that can be directly linked to the continuum stress. Let 
be the force exerted on the ith atom, and we assume that can be decomposed as 
follows. 



To our knowledge, all the empirical potentials for crystalline solids admit natural 
decompositions in this form. 

Next, within the support of (pn denoted by l^n, we have. 




(4.6) 



F„= (^„(y),f(y)). 




(4.7) 



K(y),f(y)) 



(4.8) 
(4.9) 





(4.10) 




(4.11) 




(4.12) 



^i^^n(yi)Vo + X] X] fij [Pn{yi) -^\{yj - ' "^^niYi) • (4.13) 




Here we have defined a local stress. 




(4.14) 
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This has been motivated by the expression of virial stress [80, 81]. 

The last term, (p^ (yi ) + y {jj —Ji ) • V(/?n (ji ) , can be regarded as an extrapolation 
of the basis function to the boundary of where the function vanishes. Therefore 
this term will be neglected. 

With these approximations, the sum is reduced to, 

(My), f (y)) « - E '^^'^'fin{y^)Vo. (4.15) 

In the implementation, this will be approximated by an integral, and further 
approximated by mid-point quadrature formula in each triangle, 

V CT^\/^n{y^)Vo ^ [ a {y)V ^ n{y) dj . (4.16) 

At this point, one might notice that with these approximations, the procedure in 
the CG region is the same as solving the elasto-dynamics equation, 

pouu = V-a(Vu), (4.17) 

with the stress-strain relation given by the Cauchy-Born rule [32]. 

5. Another Numerical example. Here we consider a test problem for a dis- 
location dipole in aluminum. The atomic potential used here is the EAM potential 
[31] with lattice constant 4.032A. The system under consideration is a 31) rectangular 
sample, with the three orthogonal axes being along the (110), (001) and (110) direc- 
tions respectively. The system is periodic in the third direction with period equal to 
8 times the atomic spacing. Therefore, we will treat it as a two-dimensional problem 
in the coarse-graining procedure. 

The dimension of the entire system is 45.62 nm x 48.38 nm, and the full sys- 
tem contains 384,000 atoms. Two dislocations are introduced with core position at 
(±30,0)Aand opposite Burgers vectors, given by b = ±y[110]. The initial position 
of the atoms, denoted by z^, is computed from superimposed analytical solutions of 
the anisotropic elasticity problems around edge dislocations [68, 71]. Meanwhile, a 
pure shear strain with 621 = 0.02 is applied at the remote boundary. This shear 
strain is also added to the atomic positions, and the new position will be set as the 
reference position of the atoms. Namely, we set = -\-Ezi with E being the strain 
tensor. The initial velocity is zero. For the time integration, we choose the step size 
At = 0.052880 (pico-second). 

In figure (5.1), one can see the projected position of the atoms as well as the finite 
element mesh around the two dislocation cores. In addition, in the bottom figure, we 
have indicated the selected nodes, where the subspace is extended. These nodes will 
be regarded as repeated nodes since additional degrees of freedom are introduced at 
those points. 

We first show the results from the full molecular dynamics model in Fig. 5.2. 
Since the initial position is determined from the corresponding continuum model, we 
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Fig. 5.1. Finite element nodes. Top: All the atoms in the system; Middle: the selected model 
point for a conventional Galerkin method; Bottom: Highlighting the repeated nodes at the interface. 



expect that away from the dislocation core, the atoms are almost at equilibrium at 
the beginning. Near the dislocation cores, the analytical solutions are not accurate 
due to the large deformation and loss of symmetry. This is confirmed by the numer- 
ical results: The change of displacement starts from the dislocation cores and then 
propagates into the surrounding region. 

In Fig. 5.3, we show the error of the Galerkin methods. Clearly, the conventional 
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Fig. 5.2. The snapshots of the displacement ui computed from the full molecular dynamics 
model. 



Galerkin method exhibits large error inside of the atomistic region, which is a result 
of reflection from the interface. In contrast, the extend Galerkin method introduces 
much less reflection, and offers better accuracy. The accuracy improves for further 
extensions. 
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Fig. 5.3. Snapshots of the numerical error. From left to right, T = 2,8,16, and 24. From 
top to bottom: the error from the conventional Galerkin method, the extended Galerkin method with 
L = 4, and the extended Galerkin method with L = 8 for the Krylov subspaces (4.1). 
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6. Conclusion and discussions. We have presented a new Galerkin formula- 
tion to coarse-grain a molecular dynamics model for crystalline solids. Our method 
has been designed to address the following issues: 

1. The artificial reflections at the interface. 

2. The difficulty in pre-computing the memory functions in the generalized 
Langevin equations. 

3. The stability of CG models. 

The method is referred to as an extended Galerkin method to distinguish it from the 
conventional Galerkin method and to emphasize the variational formulation. 

It is possible to extend the current formulation to finite temperature. The key is 
to satisfy the second fiuctuation-dissipation theorem [44]. One needs to introduce a 
Gaussian noise to be consistent with the approximate kernel function. This is work 
underway. 
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